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We present results from an ongoing study of mass splittings of the lowest lying states in the 
charmonium system. We use clover valence charm quarks in the Fermilab interpretation, an 
CNI improved staggered (asqtad) action for sea quarks, and the one-loop, tadpole-improved gauge 

action for gluons. This study includes five lattice spacings, 0.15, 0.12, 0.09, 0.06, and 0.045 fm, 
with two sets of degenerate up- and down-quark masses for most spacings. We use an enlarged set 
of interpolation operators and a variational analysis that permits study of various low-lying excited 
c3 states. The masses of the sea quarks and charm valence quark are adjusted to their physical values. 

This large set of gauge configurations allows us to extrapolate results to the continuum physical 
point and test the methodology. 
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Table 1: Parameters of the ensembles used in this study. Simulation light and heavy bare sea quark masses 
are denoted ;«/ and m/,. Mass-independent, tuned physical bare masses are denoted m for degenerate up and 
down quarks and m s for strange. Also shown are the approximate lattice spacing, the lattice size, and the 
number of source times used (typically four per gauge configuration). 
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1. Objectives 

The wealth of excited charmonium states discovered at the B factories presents a challenge 
for interpretation (and, in some cases, confirmation). Some states could be spin-exotic hybrids 
and some "molecular" states. In principle, lattice QCD should provide a reliable guide to the 
interpretation of these states [1, 2, 3, 4, 5, 6]. For levels above the open charm threshold, the 
treatment of multihadronic scattering states presents a technical challenge. Other challenges are 
controlling heavy-quark discretization errors and, for some levels, including annihilation effects. 
Matching theory to experiment requires a relatively high degree of precision. 

In the present study our more limited objective is to lay the foundation for future work by car- 
rying out a high-precision study of the splittings of the low-lying states. Here we use clover charm 
quarks in the Fermilab interpretation [7]. To the extent we can reproduce the known splittings, we 
test the methodology. Validation of the method also gives confidence in other studies that use the 
same fermion formulations. 

This work expands and extends our previous study [1] with clover (Fermilab) quarks and 2+1 
flavors of asqtad sea quarks. We use a large variational basis of interpolating operators, and we 
extrapolate to zero lattice spacing and physical sea-quark masses. 

2. Methodology 

2.1 Gauge-field ensembles 

We work with a large set of gauge-field ensembles generated in the presence of 2+1 flavors of 
asqtad sea quarks and a 1-loop tadpole improved gauge field [8, 9]. Parameters are listed in Table 1. 
These tuned values for the light quark masses are determined from chiral fits to masses and decay 
constants of the light pseudoscalar mesons [9]. 
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Table 2: Examples of charmonium interpolating operators used in this study. Here operators for the Tf c 
irrep are shown. In the notation below, V, generates a discrete covariant difference in direction i, Oj. = 

\etjk\ViVj, tmdM i = e ijk V i V j . 
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2.2 Interpolating operators 

As mentioned above, we use clover charm quarks with the Fermilab interpretation. Inter- 
polating operators are classified according to their cubic group irreps and their P and C quantum 
numbers. We use a large basis of operators following Liao and Manke [10] and the JLab group [11]. 
Our operators are constructed from stochastic wall sources. Averaging over stochastic sources re- 
sults in both local and smeared bilinears of the form 

ffi(x) = q{x)O iq {x) e\ (x) = q(x)OiSq{x) , (2. 1) 

where 0, is one of several operators, as illustrated in Table 2, and 5 represents covariant Laplacian 
smearing. (Only one width of smearing is included.) We do not include, however, any explicit open 
charm or charmonium/light-meson states in the list. 

2.3 Variational determination of energy levels 

In each channel (defined by the cubic group irrep and P and C), we use the standard variational 
methodology [12, 13] for determining the lowest lying states. For interpolating operators 0;, we 
define the correlation matrix 

Qjit) = (Oi(t)Oj(p)) . (2.2) 

The goal is to determine, to good approximation, the energies E n in the spectral decomposition 

v exp(-ff n Q 
C U K t )=h ZinZ jn ^ ■ ( 2 - 3 ) 

Ideally, the eigenvalues of the transfer matrix r f ~'° from time to to time t are A„ (? , to ) = exp [—E n (t — 
to)]. With a finite set of interpolating operators and a truncated list of energies, the eigenvalues 
receive contributions from higher states, which are modeled with [11] 

K(t,to) = flo exp [-E n (t - to)] +£?i exv[-E' n (t-t )] + .... (2.4) 
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The lowest energy E n for each n becomes an estimate of the nth excited state in that channel. 

2.4 Scale and charm quark mass 

The lattice scale is based on the Sommer r\ parameter. It is determined in two steps. First, we 
have detailed measurements of its value in lattice units, r\/a, from the heavy-quark potential. Then, 
r\ itself is determined on the same lattice ensembles from a mass-independent, partially quenched, 
staggered #PT analysis of the light pseudoscalar masses and decay constants. That analysis used 
fn, M n , and Mk to determine the physical bare light quark masses m and m s and r\ = 0.3 1 17(22) fm 
[14]. 

Once the lattice scale is known, we tune the charm quark mass by requiring that the kinetic 
mass M2 of the D s meson, defined through the dispersion relation, 

£ Di (p)=M 1 +p7(2M 2 ) + ..., (2.5) 

matches the experimental value. In our implementation of the Fermilab method, we do not tune 
the value of M\, the rest mass, however. Thus, the rest mass of the quark suffers from substan- 
tial discretization artifacts (except when am c <C 1). This contribution cancels, however, in the 
difference of two hadron rest masses. For that reason we report only level splittings below, e.g, 
Mi (25) — Mi(l»S). For further discussion of this point, see Refs. [1, 7]. 

2.5 Sea-quark mass effects and continuum fit model 

Although our gauge-field ensembles were generated in the presence of 2 + 1 flavors of light 
sea quarks with fixed ratios of the bare quark masses, the initial estimate that set the simulation 
strange sea-quark mass was imprecise. In the worst case, on the 0.12 fm ensembles the strange 
sea quark was some 50% heavier than our best current estimate of its physical value in a mass- 
independent scheme (see Table 1). Thus, to extrapolate any of our measurements to the physical 
point, we include terms that model the dependence of the meson masses on sea-quark masses: 

M = M + ci(2xf+x h ) + c 2 a 2 , (2.6) 

where, in the notation of Table 1 , xt = (mi —rh)/ m s and x/, = (nth — m s ) / m s . 

3. Results 

IS hyperfine splitting The hyperfme splitting of the IS level provides a demanding test of the 
methodology [16]. The extrapolated result shown in Fig. 1 is compared with the current PDG value 
and the recent BESIII value[15]. Annihilation effects, which would decrease the splitting slightly 
[17] have not been included. The extrapolated error is ±2 MeV and appears to favor the BESIII 
value. The largest contribution to the uncertainty comes from our imperfect knowledge of r\ . In 
the determination of the hyperfine splitting, this uncertainty enters twice, first in setting the charm 
quark mass, and second, in comparing the splitting with the experimental value. In this quantity 
the error is amplified, not cancelled. 

IP — IS splitting Results for the spin-averaged IP — IS splitting are shown in Fig. 1. The error 
at the physical point (including the scale error) is ±4 MeV (1%). 
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Figure 1: Left panel: 15 hyperfine splitting. Right panel: spin-averaged IP — 15 splitting. Errors on data 
points are statistical only. The blue square is the physical value including the r\ scale error. The magenta 
burst (slightly displaced) is the PDG value. The blue burst is the recent BESIII result [15]. 



Tensor Spin — orbit 
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Figure 2: Left panel: IP tensor mass combination. Right panel: IP spin-orbit mass combination 

IP tensor and spin-orbit mass combinations The left panel of Fig. 2 shows the tensor mass 
combination 

-[3M(Xci)-M(Xc2)-2M(Xco)\ 
(extrapolated error ±0.7 MeV), and the right panel, the spin-orbit mass combination 

l -[5M{Xc2)-2M{Xco)-^M{x c i)] 

(extrapolated error ±2 MeV). In a heavy-quark expansion these mass splittings arise from the 
tensor and spin-orbit terms in the effective heavy-quark potential (quark model). These results are 
sensitive to discretization errors arising from those terms. 

25 states We overestimate considerably the splitting of the 25 and 15 levels, as shown in the 
left panel of Fig. 3. We get 675 ± 6 MeV compared with the experimental value 606 MeV. The 
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Figure 3: Left panel: Spin-averaged 25— IS splittings. Note the experimental value (magenta burst) at the 
lower left. The physical open charm threshold is also shown. Right panel: y/(2S) level as a function of 
f m i„ at a sa 0.09 fm. The full energy was reconstructed by adding the measured splitting to the experimental 
spin-averaged IS level. 

same problem was seen in [1] but with less clarity. Since the 2S levels are close to the open charm 
threshold, one may speculate that by not including an explicit open charm term in the variational 
mix, we cannot get a good representation of these states [18]. In support of this hypothesis, in 
the right panel of Fig. 3, we note that the splitting of the \l/(2S) level from the IS ground state 
shows a decreasing trend as ? m in, the minimum of the fit range in A w (f), is increased. Indeed, if we 
set a higher t m [ n value for all lattice spacings and repeat the analysis, we get 651(12) MeV. Such 
behavior would result if a substantial open charm component is required but the transfer matrix 
has only a very weak mixing between closed and open charm. Such weak mixing has been known 
from string-breaking studies of the static potential [19, 20]. 

4. Conclusions and Outlook 

We can reproduce the splittings of the lowest-lying charmonium levels to a precision of a 
couple of MeV. At this level of precision, however, we fail to reproduce the 2S-1S spin-averaged 
splitting with our set of interpolators. To complete the analysis we will develop a complete error 
budget. We will next try adding explicit open charm to the variational mix. We plan, also, to study 
bottomonium. 
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